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Abstract 

We reinvestigate the dynamical behavior of a first order scalar nonlinear delay 
differential equation with piecewise linearity and identify several interesting features 
in the nature of bifurcations and chaos associated with it as a function of the delay 
time and external forcing parameters. In particular, we point out that the fixed 
point solution exhibits a stability island in the two parameter space of time delay 
and strength of nonlinearity. Significant role played by transients in attaining steady 
state solutions is pointed out. Various routes to chaos and existence of hyperchaos 
even for low values of time delay which is evidenced by multiple positive Lyapunov 
exponents are brought out. The study is extended to the case of two coupled systems, 
one with delay and the other one without delay. 
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1 Introduction 

Chaotic phenomenon has been studied extensively in several dynamical systems, including 
Chua's circuits, cellular neural networks (CNNs) and cellular neural networks with delay 
(DCNNs) during the past few decades. Dynamical systems described by delay difference- 
differential equations have many practical applications in areas such as biology, secure 
communication, economics, traffic control, neural networks, etc. For example, Mackey and 
Glass[f977] have demonstrated a model for blood production in patients with leukemia in 
terms of a delay differential equation, where the concentration of blood can vary chaoti- 
cally when request for more blood is made. Time delayed dynamical systems can also play 
a significant role in secure communications as these systems can have chaotic attractors 
with a large number of positive Lyapunov exponents when the delay time of the dynamical 
system is increased [Goedgebuer et al., 1998] and so messages can be more securely trans- 
mitted through chaotic synchronization. Time delayed dynamical systems are essentially 
infinite dimensional systems, as an infinite set of independent numbers are required to 
specify an initial condition. Such infinite dimensional systems can also be used to forecast 
the fluctuation in the share market behavior [Gwynne, 2001]. Various parameters in the 
traffic control theory are also expressed in terms of delay coordinates [Helbing, 2001]. The 
dynamics of delayed systems depends on the current state as well as on its past history 
and recursively further past history and so on. Thus CNN's with delay coordinates may 
be considered as ideal models for neural networks which act with their past knowledge. 
It is a widely accepted fact that chaos can occur in autonomous time continuous nonlin- 
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ear systems having order greater than two and in nonautonomous time continuous nonlinear 
systems with order greater than one. In time discrete systems chaos can occur even in first 
order invertible maps and in non-invertible maps with order greater than one. Recently 
Lu and He [1996] have shown that chaos can occur even in a simple scalar first order de- 
layed nonlinear (piecewise linear) dynamical system with large enough delay. Thangavel et 
al.[1998] have made preliminary study of the bifurcation scenario and controlling of chaos 
in first order and coupled DCNNs including the model of Lu and He. 

In this paper, we have revisited the first order time delayed system studied by Thangavel 
et al.[1998] with the addition of a constant external forcing term and studied its dynamics 
in a large parameter regime of the external force and time delay. In particular, we have 
identified the novel aspects of existence of a stable island in a two parameter space for 
the equilibrium solutions in the presence of time delay, the significant role of transients in 
attaining steady state solution and thereby identifying the existence of familiar routes to 
chaos and the emergence of hyperchaos even for small values of time delay in a scalar first 
order piecewise linear delay differential equation for a wide range of parameters, which 
is confirmed by the existence of multiple positive Lyapunov exponents. The study has 
also been extended to the case in which the first order scalar time delayed system is 
coupled to a second scalar system without delay. Again the existence of a stable island for 
equilibrium points is established. Existence of different bifurcation routes including type III 
intermittency route is pointed out, where we have also discussed about the transient effects. 
We have also plotted the two parameter bifurcation diagrams for both the cases to bringout 
the dynamics in two parameter space. The plan of the paper is as follows. To start with, a 
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linear stability analysis is carried out both in the presence and in the absence of time delay 
in sec. 2. In sec. 3, we have made a detailed study of the effect of transients in attaining the 
steady state solution leading to identification of typical bifurcation scenario. One and two 
parameter bifurcation diagrams for different nonlinear parameter regimes characterized by 
the piecewise linear function for different values of the parameter p 2 are studied in detail. 
We also bringout the hyperchaotic nature of the systems with the existence of multiple 
positive Lyapunov exponents. Further we have considered the addition of one more cell 
without delay in sec. 4, resulting in two cell DCNN and the results of this system are 
discussed in sec. 5. The results are summarized in sec. 6. 



2 Single cell Delayed Cellular Neural Network 

We consider the following first order delay differential equation introduced by Lu and 
He[1996], but with the addition of a constant external force as studied by Thangavel et al. 
[1998], 



< ^L = ^ ax{t ) +h f {x {t-r))+c ) (1) 
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where a and b are parameters, r is the time delay, c is a constant external force and / is 
an odd piecewise function defined as 



0. 



Mr)) = { 



x. 



0. 



x < -p 2 



■1.5x — 2, — p 2 < x < —pi 
—pi < x <p\ 



(2) 



-1.5a; + 2, Pi < x < p 2 



x > p 2 . 



Here p\ and p 2 are parameters. The form of the function f(x(r)) is sketched in Fig. ^ In 
our study we have fixed the parameter pi at p\ = 0.8 and studied the system behavior in 
a region of the external force c G [—0.15,0.15] and time delay r e [0.0,30.0] for different 
values of the parameter p 2 characterizing the piecewise linear function f(x(r)). As any 
time delayed system depends strongly on its past history, which acts with its rich past 
knowledge, this model may be viewed as a simple one cell delayed cellular neural network 
(DCNN). To start with let us consider the nature of the fixed (equilibrium) points of the 
system (1) in some detail. 



2.1 Fixed points and linear stability 

Equation (1) with the form (2) for f(x(r)) can admit upto a maximum of three fixed 
points, Xi(t) = xo,i = 1,2,3, depending upon the parameter values. A linear stability 
analysis of these fixed points under the linear perturbation x = x* + aexp(At), a < 1, for 
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Eq. (1) was carried out in [Thangavel et al., 1998]. However, the stability analysis in the 
presence of time delay was done very briefly. In this section, we bring out the existence of 
stable island in the (r, 1.56) plane for some of these equilibrium points in the presence of 
time delay in a detailed manner by following the analysis of Ramana Reddy et al. [2000] 
for the case of coupled limit cycle oscillators with time delay coupling/ time delay feedback. 
For this purpose we examine the stability nature of the fixed points for Eq. (1) both in the 
absence and in the presence of time delay r by considering the characteristic equation for 
the eigenvalue A and identifying the cases where Re(A) < for stability. We consider two 
cases. 

2.1.1 Time delay r = 

In the absence of time delay, the following fixed points can exist depending on the choice 
of parameters a, 6, c,pi and p 2 in Eqs. (l)-(2): 

a) For \x\ > P2, the fixed point is x = x$ = ^ and the characteristic equation in this region 
is A = —a. The fixed point x = x = ^ is stable for positive values of a. 

b) For — p 2 < x < —pi, the fixed point is x = xq = ^+i5b anc ^ ^he characteristic equation 
in this region becomes A = —(a + 1.56). The fixed point is stable when a > —1.5b. 

c) For \x\ < pi, the fixed point is x = xq = ^ and the characteristic equation becomes 
A = — (a + b). The above fixed point is stable when a > —b. 

d) For p\ < x < p2, the fixed point is x = xq = and the characteristic equation in 
this region becomes A = —(a + 1.56). The fixed point is stable when a > —1.56. 

Next we consider the case when time delay is present, r > 0. 
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2.1.2 Time delay r > 

a) Here again, for \x\ > P2, the fixed point remains the same as for the case r = 0, namely 
x = Xn = - and it is stable for a > 0. 

b) Next in the region, — p 2 < x < —p 1: for the fixed point x = x = ^~^ b , the characteristic 
equation becomes the transcendental equation, 



Let A = a + i(3, where a and (3 are real. Substituting this into the above equation and 
equating real and imaginary parts, we obtain equations for a and (5 as 



From the above one can easily check by squaring and adding Eqs. (4) and (5) that 



A + a + 1.56e _AT = 0. 



(3) 




(4) 



P- 1.5be- aT sm(pr) = 0. 



(5) 



= p± = ±^2.2hb 2 exp(2ar) - (a + a) 2 , 



(6) 



and 



a = —a — 



tan(/9r) 



(7) 
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Without loss of generality we choose + sign in the above equation, as the eigenvalues occur 
in complex conjugate pairs. In order to find the critical stability curves, we choose a — 0. 
Then we have 

[3\ a=0 = ±V2.256 2 - a 2 . (8) 

From Eq. (4), it follows that 



pr = icos" 1 ( j-j^ ) + 2n7r, (9) 



where n is any integer (0, ±1, ±2, ...). Consequently one can expect that the stability 
regions are confined between the set of two curves 



2rar + cos -1 (f£) 

7i(n, 1.56) = Kl - 5 bJ , 10a 

V ; V2.256 2 - a 2 y ' 

2nn - cos _1 (--^) 

r 2 (n, 1.56) = VL5 fc; . 10b 

1 ; V2.256 2 - a 2 1 ^ 



In Eq. (10a), n — 0, 1,2, ... and in Eq. (10b), n = 1,2, ... for the pair of curves we have 
chosen so that the curves have positive value of r. To start with we note from Eq. (3) that 
when the time delay r = 0, A = —a — 1.56 and so a < 0. In order to identify those curves 
for r > which encompass the stable regions, the critical curves should be the ones on 
which g > 0. From Eq. (3), 



d\ 1.56A exp(— Ar) 

= 1 - rl.56exp(-Ar) ( ' 
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and hence 



da 
dr 



a=0 



„ 1.56(z/3)exp(-i/3r) 

-Tt6 

1 — rl.5&exp(— i(3r) 
1.56/3 sin (/3r) 
5 



(12) 



where 

D = [1 - 1.56rcos(/3r)] 2 + [1.56r cos(/3r)] 2 . 

Therefore 

— — > on both n and r 2 . (13) 
ar 

The above condition implies that there can be only one stability region between r = 
line (where a < 0) and the critical curve Ti(0, 1.56) which is the closest to the line r = 0. 
We note that the condition (13) prohibits the existence of any other stable region (that is 
multistability regions) as in the case of time delayed limit cycle oscillators [Ramana Reddy 
et al., 1998; 2000], because for a second stable region to exist one requires ^ < on any 
one of the other curves (n > 1). But this never occurs in our case. The numerical plot in 
Fig. Elof the curves Ti(n, 1.56) (solid curve for n = 0, 1, 2) and ^(n, 1.56) (broken curve for 
n = 1, 2) reveals that the region between r = and r = Ti(0, 1.56) is the only stable region 
(shaded region), where on t\ > 0, which passes from negative to positive values of a, 
whereas the other curves Tiin, 1.56) < r < 7"i(n, 1.56) for n > do not satisfy the required 
stability condition and hence they are all associated with unstable regions. 
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c) As in the previous region, for the case —p\ < x < pi and for the fixed point xq = -^5, 
the characteristic equation in this region has the form 



A + a + be 



-At 



0. 



(14) 



As before, we obtain a set of critical curves as 



n(n,b) 



T 2 (n,b) 



2mr + cos- 1 (f) 



2 — n 2 



Vb' 

2wr-cos- 1 (f) 



(15a) 
(15b) 



where n = 0, +1, +2, ... in Eq. (15a) and n = +1, +2, ... in Eq. (15b). Also we get 



da 



2 n-1 



P 2 D 



a=0 



where 



D = [1 — 1.5&T cos(/?r)] 2 + [1.5&T cos((3r)f 
(3 = 6sin(/3r) 



with the same argument as for the above case, we find that there is only one stable island 
between the curves r = and ri(0, b). 

d)For the case pi < x < p 2 and for the fixed point x = x = ^tith ' ^ ne characteristic 
equation turns out to be identical to Eq. (3) and the associated critical curves also have 
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similar form as that of Eq. (10). 



3 Numerical study of the single cell DCNN 

In this section, we will present a discussion of the dynamics of the delayed cellular neural 
network (Eq. (1)) in the pseudospace (x(t),x(t — r)). Further we will also discuss the sig- 
nificant role played by transients in attaining the steady state behavior, the computational 
efficiency required for achieving such steady state solution and the nature of bifurcation 
diagrams for both low and high transients with the external forcing c as the control pa- 
rameter. We will also point out the existence of a stable island for equilibrium points in 
the two parameter (now in the time delay r - the external forcing c plane) bifurcation 
diagrams for various ranges of control parameters and nonlinearity characterized by the 
function f(x(r)) for three different values of P2 in Eq. (2). In addition, we have also cal- 
culated the Lyapunov exponents associated with the system and show that there exists 
a large parameter range over which hyperchaos exists (corresponding to multiple positive 
Lyapunov exponents). From the two parameter bifurcation diagrams, it becomes evident 
that the sizes of the hyperchaotic regime and the stable fixed point regime increase with 
time delay. 

3.1 Dynamics in the pseudospace 

The dynamics of the DCNN can be studied in a suitable phase space by plotting the numer- 
ical solution of Eq. (1) appropriately. Now, to calculate x(t) from Eq. (1) for times greater 
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than t, the function x(t) over the interval (t — r,t) must be given. Hence, for a prescribed 
continuous function x(t) on (— r, 0) one can integrate Eq. (1) using numerical methods as 
in the case of ordinary differential equations. We have numerically integrated Eq. (1) using 
Runge-Kutta fourth order integration routine with the parameters a = 1.0, b = 1.2 and 
Pi = 0.8 for three different values of P2 = 1.0, 1.33, 1.66 with the initial condition x(t) = 0.9 
in the interval (— r, 0). As the system is of first order in nature, the dynamics can be viewed 
in a pseudospace by plotting x(t) against x(t — r). (Note that the choice of x(t — r) as 
the second phase variable is arbitrary; x(t — t') can be equivalently used, where t' is an 
arbitrary time delay). One encounters typical scenario of bifurcations leading to chaos, but 
with an important difference: Transients play a crucial role and it takes a very long time 
before transients die down to attain steady state solutions. In particular, period doubling 
and inverse period doubling phenomena, besides other bifurcations, are often encountered. 
In many cases the system exhibits hyperchaotic behavior characterized by multiple positive 
Lyapunov exponents, which we have discussed below in sec. 3.4. Fig. |3] shows a typical 
pseudochaotic attractor for r = 5.0, c = 0.001 and pi = 1.33 with the initial condition as 
given above in the phase space of x(t) against x(t — 5.0). The positive maximal Lyapunov 
exponent of this attractor for the above mentioned parameter values is X ma x — 0.05461. 
It may be noted that previously it was reported by Thangavel et al. [1998] and Lu and 
He [1996] that the system for a different set of parameters exhibits only chaos (and not 
hyperchaos) for large values of time delays , that is r > 20. On the otherhand we have 
identified the parameter regimes where the system exhibits chaos and even hyperchaos for 
small values of time delay around r = 5.0. In the following we present the details. 
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3.2 Transients 

Usually dynamical systems described by ordinary differential equations can attain the 
steady state within a few thousand transients for optimal value of time step At, provided 
due care is given for numerical accuracy and suitable numerical algorithm is chosen. How- 
ever for the delay differential equation considered in this paper, such low transients (of the 
order 10 4 and more) do not seem to lead the system to their steady state solutions. In 
such cases one has to leave out more number of transients so that one can have the possi- 
bility of obtaining a clear picture of the bifurcation diagram with the existing routes. In 
fact, the effect of transients has been pointed out by Becker and Dorfler[1989] even for the 
logistic map. As for as nonlinear maps and ordinary differential equations are concerned, 
with some effort it is possible to realize whether the system has reached the steady state 
or not by constructing the appropriate bifurcation diagrams and identifying the various 
bifurcation routes by leaving sufficient number of transients before starting the analysis for 
steady state solutions. The same considerations hold good for dynamical systems modeled 
by finite number of coupled nonlinear maps and odes without time delay. 

On the other hand, when one starts studying nonlinear dynamical systems with time 
delay, effectively one considers infinite number of coupled odes which in a numerical sense 
corresponds to several hundred coupled nonlinear odes. For example, for Eq. (1), when the 
time delay r = 25.0, in the actual numerical analysis, we typically take the optimal time 
step as At = 0.05, so that we are actually solving 500 coupled nonlinear odes. However, the 
optimal time step At to be fixed depends explicitly on the nature of the system and on the 
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time delay introduced in the system. In contrast to the nonlinear odes, where either a very 
small or a relatively large time step At results in an exponentially increasing numerical 
error, for the time delayed system (1), we have found that the optimal time step At varies 
between 0.05 to 0.0001, depending on which the actual number of coupled nonlinear odes 
increase, and the number of transients to be left out is found to lie anywhere between 
2.0 — 3.0 x 10 6 . For small values of time delay r, the optimal time step At should be as 
small as possible to identify a typical bifurcation scenario. However, for large values of 
time delay r, a very small value of time step At can again be problematic as one has to deal 
with a large number of coupled odes and the time required for solving such huge number 
of coupled odes becomes huge. Also a large value of time step At does not seem to lead 
to the steady state solution, as such time steps lead to propagation of numerical errors. 
However, an appropriate optimal choice of time step At can lead to a typical bifurcation 
scenario. 

Thus the numerical analysis of time delayed nonlinear systems requires an appropriate 
time step At to be chosen, depending on which the actual number of coupled nonlinear 
odes increases, which in turn necessarily increases the number of transients to be left 
out in order to obtain the steady state solution. In such cases computing with a standard 
personal computer like a 32bit Pentium IV processor is a difficult and time consuming task. 
However, it turns out that with a high speed computers like an Alpha workstation, such a 
task becomes feasible. We have carried out our analysis with a 64 bit Alpha workstation, 
in which it took 50-60 days to plot each of the two parameter bifurcation diagrams leaving 
out sufficient transients using Runge-Kutta fourth order integration routine. 

14 



In our numerical analysis we have chosen the optimal time step as At = 0.05 for the 
value of time delay r = 25.0. The effect of transients in reaching the steady state behavior 
can be realized from the bifurcation diagrams we have obtained for the parameter values 
a = 0.16,6 = 0.2 and for two different values of p%. Fig. |3k shows the one parameter 
bifurcation diagram for the above mentioned parameter values with pi = 1.33, when we 
leave transients of the order 1.0 — 2.0 x 10 4 , which is sufficiently large in the case of 
nonlinear maps and odes, leaving the impression that the steady state has been reached 
with a bifurcation scenario that is quite complicated and atypical. In contrast, Fig. Hb 
shows the bifurcation diagram for the same value of the parameters except that now the 
transients left out are of the order of 1.0 x 10 5 , which shows a typical bifurcation scenario. 
Similarly Fig. shows the bifurcation diagram for the same values of parameters as above 
for r, a and b with p 2 = 1.66 for the transients of the order 1.0 x 10 5 , which shows a very 
complex structure, whereas Fig. shows the bifurcation diagram for the transients of 
the order 1.4 x 10 6 , in which case still there exists some stray points near the bifurcation 
regions, which indicates that the system requires still more transients to settle to its steady 
state. The complexity increases with time delay and hence the number of transients also 
increases, which in turn increases the computing time enormously. This is also evident 
from Fig- El where the maximum value of x(t), XMAX, is plotted against the time delay r, 
for the parameter values a = 1.0, b = 1.2, c = 0.001 and p 2 = 1.33. Fig. shows that even 
for transients of the order 1.0 x 10 5 , the bifurcation diagram is still in an unsettled form for 
larger values of time delay r. On the other hand in Fig. the number of transients is of 
the order 2.5 x 10 5 , wherein a clear bifurcation scenario has emerged. We have verified the 
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role of transients to attain the steady state solutions for various values of time step At and 
time delay r for few other delay differential equations also, including those discussed in the 
work of Ramana Reddy et al., [2002]. The results of our detailed numerical investigation 
for Eq. (1) are tabulated in Table I, where we have indicated the number of transients to 
be left out for given values of At in order to attain the steady state solutions. 



Table I: Effect of step-size on transients 



Value of r 


Optimal value of At 


Number of coupled 
differential equations 


Number of transients 
to be left out 


5 


0.002 


250 


3.0-5.0 x 10 5 




0.001 


500 


5.0-6.0 x 10 5 


10 


0.002 


5000 


1.0-2.0 x 10 6 




0.001 


10000 


> 2.0 x 10 6 


25 


0.05 


500 


1.0 x 10 5 - 2.0 x 10 6 




0.005 


5000 


1.0-2.0 x 10 6 




0.0025 


10000 


> 2.0 x 10 6 



As the time delay increases, the size of the delay loop (that is to be updated on every 
iteration) required to maintain the delay variable in numerical simulation also increases, 
which in turn increases the computational effort. Another important point to note in 
connection with the effect of transients is the following. In the above we have discussed 
the effect with reference to a scalar delay differential equation only. On the other hand, in 
recent times, there has been considerable interest to study the dynamics of multiply coupled 
delayed neural networks in biological systems [Campbell et al., 1999; Earl & Strogatz, 2003]. 
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In order to examine the actual steady state behavior of such networks, which effectively 
corresponds to several thousands of coupled nonlinear odes, the time required to obtain 
steady state solution increases enormously as each one of them have separate time delay 
and intrinsic transient behavior and it becomes crucial to study the effect of the later. 

3.3 One and two parameter bifurcation diagrams 

All the bifurcation diagrams shown in this paper are plotted after leaving out a very large 
number transients of of the order of 1.2 — 3.2 x 10 6 . Fig. EJo shows the one parameter 
(c) bifurcation diagram characterized by the piecewise linear function f(x) in Eq. (2) with 
Pi = 0.8, P2 = 4/3, t = 25.0 for c G [—0.1,-0.05]. It clearly exhibits period doubling 
scenario; however there exists a sudden distortion around the value of c = —0.0785, the 
cause of which remains unexplained by any standard type of bifurcation. (We wish to 
point out that the numerical analysis in [Thangavel et al., 1998] actually corresponds to 
the value of p 2 = 1.0 and not p 2 = 4/3 as stated in that paper). Similarly Fig. shows 
the bifurcation diagram with the values pi = 0.8, p 2 = 5/3, r = 25.0 within the range 
of c G [—0.15, —0.136], exhibiting reverse period doubling route to chaos in the range of 
c G [—0.136,-0.144] that includes a clear band merging crises and anti-monotonicity in 
the range of c G [-0.15, -0.144]. 

The two parameter bifurcation diagram for r G [0,30] and c G [—0.16,0.16] when 
p 2 = 1.0 is shown in Fig. [7k., which shows the behavior of the DCNN for the combined 
phase space of parameters r and c. The following colour codes are used to represent 
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various regions, period-1 region -red, period-2 region - green, 3-blue, 4-yellow, 5-magenta, 
6-cyan, 7-gray, 8-copper, chaos-black and the fixed points-white. The white region in the 
two parameter bifurcation diagram corresponds to stable regions. Fig. [TJd shows the two 
parameter bifurcation diagram for the same parameters as in Fig. [7b,, except that now 
P2 = 4/3. We use the same colour codes as in Fig. UJi for all the two parameter bifurcation 
diagrams in this paper. Similarly Fig. [7b shows the two parameter bifurcation diagram for 
the same range of the control parameters and the same values of the parameters a and b 
in Eq. (1) except that now p 2 = 5/3. By comparing the Figs. [7k,, U3p and [7b, we can see 
that the stable fixed point regions and the chaotic regions increase for small changes in 
the nonlinear parameters p\ and p 2 - Further we infer from these figures that as the delay 
increases the chaotic nature of the system also increases and thereby contributing to the 
hyperchaotic nature of the system characterized by multiple positive Lyapunov exponents. 

3.4 Lyapunov exponents and hyperchaotic regimes 

One of the interesting aspects of the dynamics associated with Eq. (1) is the existence 
of hyperchaos in a single first order scalar equation with time delay even for small values 
of the latter for suitable values of other system parameters, whereas mention was made 
in the earlier papers [Lu & He, 1996; Thangavel et al., 1998] that the system Eq. (1) 
exhibits only simple chaos even for large values of time delay. As Eq. (1) is a first order 
delay differential equation, the usual procedure for calculating Lyapunov exponents is 
not applicable. However, the simple idea of approximating a single scalar differential 
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equation with delay by an N- dimensional discrete mapping [Farmer; 1982], which is of 
infinite dimensions, facilitates one to calculate the Lyapunov exponents. To stimulate the 
behavior of such systems on a computer it is necessary to approximate the continuous 
evolution of an infinite dimensional system by a finite number of elements whose values 
change at discrete time steps. In this manner a continuous infinite dimensional dynamical 
system is replaced by a finite dimensional iterated mapping. This method was proposed 
originally by Farmer [1982] to calculate the Lyapunov exponents for delay systems. As the 
delay parameter is increased, for most parameter values the dimension increases and the 
attractor generally becomes more complicated, thereby contributing to the hyperchaotic 
nature of the system, which gets confirmed by the increasing number of positive Lyapunov 
exponents. The first ten maximal Lyapunov exponents, for the parameter values a = 
1.0, b = 1.2, c = 0.001, pi = 0.8, p 2 = 1.33 and r G (2,29) is shown in Fig. H where 
it is evident that the number of positive Lyapunov exponents increases with time delay 
r. Almost the entire black regime in the two parameter bifurcation diagrams Fig. [7| for 
large values of r is characterized by multiple positive Lyapunov exponents contributing to 
hyperchaotic nature of the system. 
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4 Stability analysis and chaotic dynamics of two cell 



DCNN 

We now add one new cell without delay to the Eq. (1) to obtain the following coupled 
equation 



dX ^ = -ax(t)+bf(x(t-r)) + dy(t), (16a) 



dt 
dy(t) 
dt 



= -cy(t) + ex{t), (16b) 



where c and e are additional parameters and the function f(x) and other parameters are 
the same as already defined in Eq. (1). As in the case of single cell DCNN, we have 
fixed the parameter p\ at pi = 0.8 and studied the system behavior in a range of variable 
parameter c G [0.0, 1.4] and time delay r G [0.0, 30.0] for different values of the parameter 
p 2 characterizing the piecewise linear function f(x). Eq. (16) can also be viewed as a two 
cell DCNN, one cell acting with its past knowledge coupled with another cell which does 
not have any past knowledge. Now let us consider the nature of the fixed points of system 
(16). 



4.1 Fixed points and Linear stability 

We will bring out the existence of stable island in two cell DCNN by linear stability analysis 
in the presence of time delay r, as in the case of single cell DCNN. For this purpose, we 
examine the stability nature of the fixed points of Eq. (16) both in the absence and in the 
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presence of time delay r. 



4.1.1 Time delay r = 

In the absence of time delay, the following fixed points can exist depending on the choice 
of parameters a, 6, c, d, e,pi and p 2 in Eqs. (16) and (2). 

a) For \x\ > p 2 , the fixed point is (x,y) = (x ,y ) = (0, 0) and the characteristic equation 
in this region is A 2 + (a + c)A + ac — de = 0. The fixed point (x, y) = (0, 0) is stable when 
a + c > and ac > de. 

b) For -p 2 < x < - Pu the fixed point is (x,y) = ( de _ (a + c L5fe)c , de _ (a 2 + L5fe)c ) and the charac- 
teristic equation in this region becomes A 2 + (c + a + 1.56)A + (ca + 1.56c — de) = 0. The 
fixed point is stable when c + a + 1.56 > and ac + 1.56c > de. 

c) For |x| < pi, the fixed point is (x,y) = (0,0) and the characteristic equation becomes 
A 2 + (c + a — 6) A + (ca — be — de) =0. The above fixed point is stable when c + a — 6 > 
and ac > be + de. 

d) For pi < x < pa, the fixed point is (x,y) = ( de _^l bb)c , de -[ a +i. 5b)c ) and the character- 
istic equation in this region becomes A 2 + (c + a + 1.56) A + (ca + 1.56c — de) = 0. The fixed 
point is stable when c + a + 1.56 > and ac + 1.56c > de. 

Next we consider the case when time delay is present, r > 0. 

4.1.2 Time delay r > 

a) For \x\ > p 2 , the fixed point and its stability remains the same as for the case r = 0.0 

b) Next in the region -p 2 < x < -p u for the fixed point (x,y) = ( de _ (a + c L5fe)e , de _ (a + e 1 . 5b)c ) 
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the characteristic equation is 



(c + A) (A + a + 1.56exp(-Ar)) - de = 0, 



(17) 



which is a transcendental equation with infinite number of solutions. Let A = a+i(3, where 
a and j3 are real. Substituting this in the above equation and equating real and imaginary 
parts, we obtain two equations as 

(c + a)(a + a + 1.56 exp(— or) cos(/3r)) — de — /3 2 + 1.5/36 exp(— ctr) sin(/3r) = 0, (18a) 
(c + a)(/3 - 1.56exp(-ar)sin(/3r)) + j3{a + a + 1.56exp(-ar) cos(/3r)) = (18b) 

In order to find the critical stability curves, we choose a = 0. Then we have 



Multiplying the above two equations ()19aj) and (|19bj) with c and /3, respectively, adding 
we obtain 



1.56c cos(/3r) + ca - de - /? 2 + 1.5/36 sra(/3r) = 0, 



(19a) 



1.56/3 cos(/3r) + c/3 + /3a - 1.56c sin(/3r) = 0. 



(19b) 



1.56 cos 03r)(c 2 + /3 2 ) + a(c 2 + /3 2 ) - cde = 0. 



(20) 
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Now multiplying the equations ()19a|) and (j!9b|) with (3 and c, respectively, subtracting the 
resulting equations, we obtain 



1.5bsm(/3r)(c 2 + f3 2 ) - (3(c 2 + (3 2 ) - f3de = 0. (21) 

Squaring and adding the equations (J20|) and (|21|). rearranging them, we obtain the following 
cubic equation for (3 2 , 

X 3 + uX 2 + vX + w = 0,X = /3 2 , (22) 
where the constants are given as 

u = a 2 +2(de+c 2 )-2.25b 2 ;v = (de+c 2 ) 2 -2ac(de-ac)-A.5b 2 c 2 ;w = c 2 (de-ac) 2 -2.25b 2 c 4 . 



From Eq. (|2T)|) . we obtain 



i / cde — a(c 2 + 3 2 ) , 



where n is any integer (0, ±1, ±2, ...) and the value of (3 can be obtained by solving the 
Eq. (}2*2*|) for (3 2 . Consequently the stability regions are confined between the set of curves, 



2mr + cos- 1 t cde - a{c2+ P 2) 



, 1.5b(c 2 +/3 2 ) , 

-,.("■ l.->h) = j / (2 La) 

2/^ -cos" 1 



, 1.56(e 2 +/3 2 ) . 

r 2 (n, 1.56) = j 7 , (2 lb) 
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where n — 0, +1, +2, ... in Eq. (j24aj) and n = +1, +2, ... in Eq. (|24b|) . In order to check 
whether the region enclosed by the curve T\ and Ti for every value of n forms stable islands 
in (r, 0) plane, one has to examine the sign of on T\ and r 2 as we have done for the case 
of single cell DCNN. We have found that there exists only one stable island between r = 
and r = Ti(0,/3) curve (shaded region) as shown in Fig. El 

c) In the region \x\ < p±, for the fixed point (x,y) = (0,0), the characteristic equation is 



(c + A) (A + a + 6exp(-Ar)) - de = 0, 



is a transcendental equation with infinite number of solutions. Proceeding in the same way 
as for the case — P2 < x < — p%, we obtain a set of critical curves as 



2mr + cos 1 ( ] 



ti (n, b) = p (25a) 



2wk — COS 



_1 / cde-a(c 2 +l3 2 ) \ 



r 3 (n, b) = ^— (25b) 



where n = 0, +1, +2, ... in Eq. (j25aj) and n = +1, +2, ... in Eq. (|25b|) . As for the previous 
case we found that there is only one stable regions between r = and r = ti(0, (3). 
d)For the case p 1 < x < p 2 and for the fixed point (x,y) = ( de _[ a +t 5b)c , de -[ a +t 5b)c ) > the 
characteristic equation turns out to be indentical to Eq. (17) and the associated critical 
curves also have similar form as that of Eqs. ()24a|) and (|24bj) . 
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5 Numerical study the two cell DCNN 

In this section, we will discuss the dynamics of the two cell DCNN defined by Eq. (16) 
through numerical analysis and the transient effects involved in it by studying the bifur- 
cation diagrams for both low and high transients with c as control parameter. The one 
parameter bifurcation diagram in one of the regimes has the signature of type III inter- 
mittent behavior. We have also found that the dynamics at this transition possesses type 
III intermittent characteristic scaling behavior. We have also plotted the two parameter 
bifurcation diagrams in the (r, c) plane for three different values of P2, from which also we 
have pointed out the existence of stable island for equilibrium points for suitable choice of 
other parameter values. 

5.1 Transients 

As we have already discussed much about the role of transients in the case of single cell 
DCNN, in Sec. 3.2, it is not necessary to discuss it once again in detail for the present 
case. So we point out only the number of transients required to obtain the bifurcation 
diagrams which we have shown in Figs. El and ^2 for the system(16), which enables 
one to realize the effect of transients. As discussed in the case of single cell DCNN, the 
transients predominates evolution of the system(16) to its steady state solutions. Fig. ITTik 
shows the one parameteriii bifurcation diagram for transients of the order 2.0 x 10 5 for the 
parameter value P2 = 1.33, whereas Fig. ITUb shows the one parameter bifurcation diagram 
for transients of order 1.75 x 10 6 for the same parameter values. The stray points in the 
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neighbourhood of bifurcation points are due to the transient effects, which suggests that 
it requires still more iterations need to be left out. Fig. fTTk shows the one parameter 
bifurcation diagram for transients of the order 2.0 x 10 5 for the parameter value p 2 = 1.33, 
whereas Fig. ITTb shows the bifurcation diagram for transients of the order 2.0 x 10 6 for 
the fixed values of other parameters. 

5.2 One and two parameter bifurcation diagrams 

We have integrated Eq. (16) with the parameters a = 0.16, b = 0.2, d = 0.2, e = 0.2 and c 
as variable parameter with the initial conditions x(t) = 0.9 for t G (—25, 0) and y = 0.8. 
Fig. ITUb shows the one parameter bifurcation diagram in the nonlinear parameter regime 
characterized by the function f(x(r)) with p\ = 0.8, p 2 = 4/3 and r = 25.0, which shows 
the period-3 doubling bifurcation route to chaos. Fig. ITTb shows reverse period doubling 
in the range c G [0.32,0.345], and at the critical value of c cr u = 0.345969 the system 
exhibits intermittent transition to chaos followed by reverse period 5 doubling in the range 
c G [0.346, 0.365] interspersed by periodic windows for the same parameter values as above 
except that now p 2 = 5/3. At the intermittent transition, the amplitude variation loses 
its regularity and a burst appears in the regular phase as shown in Fig. This behavior 
repeats as time increases as observed in the usual type III intermittent scenario. The 
duration of laminar phases is random during the transition and finally results in chaotic 
oscillations by increasing the critical value of the control parameter c. The plot of mean 
laminar length < I > as a function of the parameter / = {c CT u — c) is shown in Fig.^3 where 
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c cr u is the critical value of the parameter for the occurrence of the intermittent transition. 
The phase space trajectories reveal a power law relationship of the form < / >= f~ a 
with the estimated value of a = 0.871. This analysis confirms that the trajectories at the 
critical value of c is associated with standard intermittent dynamics of type III described 
in the work of Pomeau Sz Manneville [1980] and Schuster [1988]. Figs. ITU shows the global 
bifurcation diagrams when p 2 = 1.0,4/3 and 5/3 respectively. The colour codes used 
here are the same as in the previous case of single cell DCNN. By comparing the global 
bifurcation diagrams, we can realize that the stable island (fixed point) and the chaotic 
region increases for small changes (increase) in the parameter p 2 . In addition the chaotic 
nature of the system increases with delay. 

6 Conclusion 

We have studied in detail the linear stability nature of equilibrium points, both in the 
absence and in the presence of time delay, of the DCNN systems and identified the existence 
of stable island in a single scalar piecewise linear differential equation with inherent time 
delay. The significant role of transients in attaining the steady state behavior has also been 
brought out. One and two parameter bifurcation diagrams for specific choice of parameters 
have been presented in support of our results. Hyperchaotic nature of the system, even 
for small time delay in a simple first order piecewise linear differential equation ensures 
that such systems can be used for highly efficient secure communication by constructing 
equivalent electronic circuits. Therefore based on our results the respective circuits with 
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suitable delay can be designed and can be tested for specific applications. These studies 
have also been extended to coupled systems with time delay and the generic nature of the 
results established. 
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Figure 1: The form of the piecewise linear function f(x(r)) in Eq. (1) for the values 
Pl = Q.8,p 2 = 1.33. 



Figure 2: Curves of Eq. (10a) and Eq. (10b) in the range of —p\ < x <p\. The solid curve 
represents T\ for n = 0, +1, +2 and broken curve represents t 2 for n = +1, +2. The region 
enclosed between the line r = and the curve r = t\ (0,1.56) is the only stable island 
(shaded region). 



Figure 3: Pseudochaotic attractor for r = 5.0, a = 1.0,6 = 1.2, c = 0.001, pi = 0.8, p 2 = 
1.33, x(t) = 0.9, t E (-5,0). 
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Figure 4: Bifurcation diagrams of single cell DCNN for the parameters values a = 0.16, 6 = 
0.2 and r = 25.0 when p 2 = 1.33 (a) for transients of the order 1.0 — 2.0 x 10 4 and (b) for 
transients now of the order 1.0 x 10 5 . 

Figure 5: Bifurcation diagrams of single cell DCNN for the parameter values a = 0.16, 6 = 
0.2 and r = 25.0 when p 2 = 1.66 (a) for transients of the order 1.0 x 10 5 and (b) for 
transients of order 1.4 x 10 6 . 

Figure 6: Bifurcation diagrams of single cell DCNN for the parameter values a = 1.0, b = 
1.2, c = 0.001, p 2 = 1.33 and r G (2, 29) (a) for transients of the order 1.0 x 10 5 and (b) for 
transients of the order 2.5 x 10 5 . 

Figure 7: Global bifurcation diagrams of single cell DCNN for r G (0.5,30) and c G 
(—0.16,0.16). The following colour codes are used to represent various regions, period-1 
region -red, period-2 region -green, 3-blue, 4-yellow, 5-magenta, 6-cyan, 7-gray, 8-copper, 
chaos-black and the fixed points-white, (a) p 2 = 1.0, (b) p 2 = 1.33 and (c) p 2 = 1.66. 

Figure 8: The first ten maximal Lyapunov exponents for the values a = 1.0,6 = 1.2, c = 
0.001, p 2 = 1.33 and r G (2,29). 

Figure 9: Curves of Eq. (|24a|) and Eq. (j24b|) in the range of —p\ < x < p\. The solid curve 
represents T\ for n = 0, +1, +2 and broken curve represents r 2 for n = +1, +2. The region 
enclosed between the line r = and the curve r = t% (0,1.56) is the only stable island 
(shaded region). 

Figure 10: Bifurcation diagrams of coupled cell defined by Eq. (16). (a) For the parameters 
values a = 0.16, 6 = 0.2, d = 0.2, e = 0.2 and r = 25.0 for p 2 = 1.33 and for transients of 
the order 2.0 x 10 5 and (b) For the same parameter values as in (a) except for the transients 
of the order 1.75 x 10 6 . 

Figure 11: Bifurcation diagrams of the two cell DCNN (a) For the parameters values 
a = 0.16, 6 = 0.2, d = 0.2, e = 0.2 and r = 25.0 when p 2 = 1.66 for transients of the order 
2.0 x 10 5 and (b) For the same parameter values as in (a) except for transients of the order 
2.0 x 10 6 . 

Figure 12: Intermittent behavior at the parameter values a = 0.16, 6 = 0.2, d = 0.2, e = 0.2 
and r = 25.0 when p 2 = 1.66 for critical value of the parameter c cr u = 0.345969. 

Figure 13: Mean laminar length < / > versus / = c cr u — c. 
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Figure 14: Global bifurcation diagrams of coupled cell DCNN for r e (0.1,30) and c G 
(0.3, 1.4). The colour codes are the same as in Fig. 7. (a) p 2 = 1.0, (b) p 2 = 1-33 and (c) 
p 2 = 1-66. 
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